Quasicrystalline Order in Binary Dipolar Systems 
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Motivated by recent experimental findings, we investigate the possible occurrence and character- 
istics of quasicrystalline order in two-dimensional mixtures of point dipoles with two sorts of dipole 
moments. Despite the fact that the dipolar interaction potential does not exhibit an intrinsic length 
scale and cannot be tuned a priori to support the formation of quasicrystalline order, we find that 
configurations with long-range quasicrystallinity yield minima in the potential energy surface of the 
many particle system. These configurations emanate from an ideal or perturbed ideal decoration 
of a binary tiling by steepest descent relaxation. Ground state energy calculations of alternative 
ordered states and parallel tempering Monte-Carlo simulations reveal that the quasicrystalline con- 
figurations do not correspond to a thermodynamically stable state. On the other hand, steepest 
descent relaxations and conventional Monte-Carlo simulations suggest that they are rather robust 
against fiuctuations. Local quasicrystalline order in the disordered equilibrium states can be strong. 



I. INTRODUCTION 

Since the surprising discovery of sharp diffraction im- 
ages with non-crystallographic symmetry in some rapidly 
quenched metal alloys by Shechtman et al. in 1984 
quasicrystalline structures have attracted lastly growing 
interest as an alternative type of structure of solid mat- 
ter. As a distinctive feature, these quasicrystals pos- 
sess long-range positional order in combination with a 
crystallographically 'forbidden' (e.g. fivefold) point group 
symmetry, which necessarily means aperiodic order. By 
now, many systems with a quasicrystalline order have 
been identified in nature, most of them being ternary 
or, in a few cases, binary alloys (for a review see 2]). 
Very recently, quasicrystalline structures with fundamen- 
tal building blocks much larger than single atoms have 
been found in micellar phases of dendrimers (tree-like 
molecules) Q • They represent a new mode of organiza- 
tion in soft matter and are interesting in connection with 
photonic bandgap materials Q and photonic quasicrystal 
lasers 0- 

On the theoretical side, the formation of quasicrystals 
could be reproduced in simulations of binary mixtures 
with hard sphere or Lennard- Jones potentials 0, ■ Iii 
these simulations the quasicrystalline structures were sta- 
bilized by tuning the distances corresponding to the min- 
ima of the interaction potentials to match the specific 
particle-particle distances. 

Recently, some evidence for local patters with fivefold 
symmetry was found in two-dimensional binary mixtures 
of superparamagnetic colloidal particles 0|- 
Hence the intriguing question arises: Can there be qua- 
sicrystalline long-range order in a binary dipolar system 
despite the fact that the dipolar interaction potential 
does not possess tunable intrinsic length scales? 
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Here, we tackle this problem by investigating the oc- 
currence of quasicrystalline order in binary mixtures of 
point dipoles with varying dipole strengths. As a two- 
dimensional reference structure with fivefold symmetry, 
we use the prominent rhombic binary tiling and deco- 
rate it with two types of dipoles. Then we let it relax 
mechanically and find that for a certain range of dipo- 
lar strength ratios D the mechanically stable configura- 
tion preserves the long-range quasicrystalline order (see 
sec. IH Cl below). The final configuration corresponds to 
a local minimum in the potential energy surface of the 
many-particle system. However, its global features are 
still undetermined. It could correspond to the ground 
state, a thermodynamically stable state within a cer- 
tain parameter regime (of temperature, mixing ratio and 
dipolar strength ratio), or a metastable state that, be- 
low some freezing temperature, becomes separated from 
other states by free energy barriers which are infinite in 
the thermodynamic limit of infinite system size. This 
is reminiscent to the behavior of mean-field spin-glass 
models jl^. In both cases, one would expect the dipo- 
lar system to evolve in time at some finite temperature 
and ultimately build up permanent long-range quasicrys- 
talline order if its initial configuration belongs to the 
attraction basin of the quasicrystalline state. On the 
other hand, the mechanically relaxed quasicrystal struc- 
ture could correspond to a metastable state from which 
the system escapes when it surmounts a finite free energy 
barrier. Nevertheless, within the last scenario, it would 
be interesting to see, whether quasicrystalline ordering 
exists locally in the disordered equilibrium state of the 
system. 

To evaluate the global features as discussed in the last 
paragraph, we perform a number of investigations. We 
first compare in sec. lllll the energy of several plausible al- 
ternative ground state structures with the energy of the 
quasicrystalline state. Then, we assess in sec. II VI the sta- 
bility of the quasicrystalline structure, for varying dipole 
strength ratio, with respect to small random perturba- 
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FIG. 1: The binary tiling decorated by strong (A) and weak 
(B) dipoles. The angles enclosed by the edges of its two rhom- 
bic building blocks are all multiples of 7r/5, which gives rise 
to the non-translation invariant 10-fold symmetry. 



ments mi and m2 at a distance r is given by 



^0 mim2 
47r 



(1) 



where /xq is the vacuum permeability. In the following 
we switch to dimensionless quantities. As unit length 
we choose the edge a of a rhomb, i.e., the distance of 
two neighboring A- and B-dipoles in the reference struc- 
ture. As unit of energy (and temperature T) we choose 
the interaction energy {hq/Att) ■ {niAmB / a?) of two such 
dipoles. After introducing the ratio D = mx/m^ of 
dipole strengths and me, the interaction potentials of 
two A-dipoles, an A- and a B-dipole, and two B-dipoles 
at distance r are 



Eaa = Dr 
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^'ab = r 



Ebb = D'^-"" . 
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tions of the particle positions in the ideal decoration. 
Finally, in sec. we perform Monte Carlo simulations 
to analyze the characteristics of quasicrystalline order at 
finite temperatures. Since for the thermalization of the 
system conventional Monte Carlo techniques turned out 
to be ineffective, we applied a parallel tempering protocol 
to reach thermodynamic equilibrium. 



II. MODEL 
A. System parameters and order parameter 

Corresponding to the experimental situation 0, [To| . 
we consider particles moving in a plane with their dipo- 
lar moments all pointing in the same direction perpen- 
dicular to the plane. In experiments this can be realized 
by letting superparamagnetic particles float on a liquid 
meniscus in an external magnetic field. 

To investigate quasicrystalline ordering, we choose as 
a reference structure the two-dimensional binary tiling, 
which consists of two types of rhombs put together by cer- 
tain matching rules |l3| . It is a typical example for a qua- 
sicrystalline pattern with fivefold symmetry (see fig.ni. 
A crystal results from a periodically repeated unit cell 
decorated by atoms. Here we obtain the ideal quasicrys- 
talline reference structure by decorating the rhombs as 
illustrated in fig. ^ and described in ref. so that it 
consists of A^A strong A and A'^b weak B dipoles. The 
mixing ratio x = Na/{Na + Nb) is t/(2 + t) 0.447, 
where t — {1 + \/5)/2 is the golden mean. A picture of 
the final structure is shown in fig. ^ The decoration was 
already used in previous simulations of binary Lennard- 
Jones and hard sphere systems 0, SH- However, in 
contrast to these simulations, the potential of the dipo- 
lar system cannot be optimized in an obvious way to 
support the formation of the quasicrystalline structure. 

The interaction energy of two parallel magnetic mo- 



The dimensionless form clarifies that the dipole strength 
ratio D is the only additional parameter in our problem 
besides the temperature T and the (fixed) mixing ratio 
X. This is a direct consequence of the scale-free dipolar 
interaction potential. 

In order to quantify the degree of quasicrystalline or- 
der, we define the order parameter 
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^(r-max - rj^k) exp(i • 10 aj.k) , (3) 



where 9{.) is the conventional step function {9{x) = 1 
for a; > and zero else), Np = J2 



the number of pairs of dipoles with distances rj^k smaller 
than Tinax = 1.15, and aj,k = <(rj_fc, e) is the bond angle 
between the pair vector Vj^k and an arbitrary but fixed 
direction e. Since the bond-angles in the binary tiling 
are all multiples of 7r/5 (cf. fig. = 1 in the ideal 
quasicrystalline structure, while (j) — for a system with- 
out tenfold bond orientational order. The value rmax is 
slightly smaller than the distance of the two A parti- 
cles in the fat rhomb shown in fig. ^ Eliminating these 
pairs from the sum in eq. (jJJ, on the one hand, allows 
us to encompass slightly displaced "nearest neighbors" 
but, on the other hand, guarantees that bond angles dif- 
ferent from multiples of 7r/5 are not counted in the ideal 
quasicrystalline configuration. 

In addition to we also consider the n-fold local 
bond orientational order parameters 



2N ^ N, \^ 



rj,k) exp(i -naj^k) , (4) 



where Nj — ^(rinax — fj.fe) is the number of neighbors 
of dipole j. Since the absolute value is taken before av- 
eraging over all N dipoles, the 4>n are sensitive to the n- 
fold symmetric arrangement of nearest neighbors around 
any dipoles, but insensitive to the spatial variation of the 
bond orientations. Naturally, we have 0io > (t>. 
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B. Rational approximants 

Due to the missing translational invariance, standard 
periodic boundary conditions cannot be imposed to ideal 
quasicrystalline systems. To resolve this problem, it is 
convenient to use rational approximants of quasicrystals 
as described in . A rational approximant is a rectan- 
gular part of the ideal quasicrystalline structure which is 
chosen such that it may be exposed to periodic bound- 
ary conditions without too much distorting the local qua- 
sicrystalline ordering. 

From a number of different rational approximants, we 
mainly used a small one with 890 dipoles (398 A, 492 B 
dipoles) and size 24.80 x 29.15, and a large one with 1700 
dipoles (760 A, 940 B dipoles) and size 34.27 x 40.29. 
For both approximants, we find the order parameter 
(f) ~ 0.9998 close to the ideal value (p — \. We have 
no indication that the type of approximant is decisive for 
the results obtained below. 

To speed up the computer simulations, we store the 
dipolar interaction energies (and the forces) in a large 
matrix by discretizing the set of possible distances r^j 
and use a linear interpolation scheme between the matrix 
entries. As an advantage, the matrix has to be calculated 
only once for all the simulations, and the computer mem- 
ory access is usually much faster than the repeated eval- 
uation of the original mathematical expressions. Details 
of the method and how to take into account the periodic 
boundary conditions are outlined in the appendix. 

As an example, for the two approximants described 
above we use a 4044 x 4754 matrix for the energies and 
a 2022 X 2377 matrix for the forces. We find the under- 
lying discretization of space for the energies or forces on 
a length scale of order 10~^-10~'^ to have no noticeable 
influence on our results. 



C. Mechanical equilibrium 

When considering the ideal quasicrystalline structure 
as a potential (meta-)stable state, the first question one 
should ask is whether this structure can be mechanically 
stable, i.e., whether the net forces on all the dipoles bal- 
ance out or, in other words, whether the structure is a 
local or global minimum of the systems' potential energy 
[l^. While in our case for common periodic lattices, this 
is obvious from simple symmetry arguments, the situa- 
tion is considerably more involved in the quasicrystalline 
structure, which is another instance of the peculiarities 
of this unusual type of ordering. 

Note that in an infinite quasicrystal two different posi- 
tions are never exactly equivalent. So the calculation of 
the long-range interaction with its surrounding dipoles 
can, strictly speaking, not be based on a finite piece of 
the ideal structure. Practically, however, as the dipolar 
interactions or Q decay as with the distance r, 
the potential felt the dipoles will largely be dominated 
by their local surroundings (see also next section). In 



this context, it is interesting to note that in the infinite 
binary tiling any arbitrarily large piece of it repeats (up 
to rotations) within a distance of the order of its size (see 
e.g. 113).^ 

In the ideal structure, the net forces acting on a dipole 
converge to a relative precision of 10^^ when taking 
into account neighboring particles up to a distance of 
Tinax — 25. The resulting force components clearly show 
that there is no dipole strength ratio D which would make 
the ideal structure a minimum of the potential energy 
surface. On the other hand, when we let the ideal qua- 
sicrystalline structure relax via the method of steepest 
descent into a local potential minimum, we find that in 
the range 4 < D < 6.5 the positions of the dipoles are 
only very slightly shifted. Based on the steepest descent 
algorithm, however, we cannot state simple systematic 
rules for the decoration of the binary tiling such that the 
dipoles assume an exact mechanical equilibrium. Due to 
this fact and since the energy per dipole of the relaxed 
configurations and the ideal quasicrystal are almost the 
same, we retain the decoration of the binary tiling as 
model reference structure. 



III. GROUND STATE CALCULATIONS 

Considering the quasicrystalline structure as a poten- 
tial equilibrium state, we next investigate how this struc- 
ture compares energetically to plausible alternative or- 
dered states. To this end, we examine the energy per 
dipole in a binary dipolar system with the mixing ratio 
X = t/(2 -f r) and the total number density of dipoles 
p ~ 1.231, which are the values of the ideal quasicrys- 
talline structure. The dipole strength ratio D then re- 
mains as the only free parameter. 

To determine the energy per dipole Ei, in the ideal 
quasicrystalline structure, we calculate the average en- 
ergy of about 10^ dipoles within a maximum distance 
''max = 300 and extrapolate for r^ax oo view 
of the discussion in the previous section, we should men- 
tion that the value of E^, may be up to a few per mille 
higher than the energy per dipole in the slightly 'relaxed' 
quasicrystal. Note that the conclusion from this section 
will not alter due to such differences. 

The irrational mixing ratio x requires alternative or- 
dered structures based on regular lattices to be phase- 
separated, i.e., to be a combination of two distinct lat- 
tices, one with a higher mixing ratio than in the qua- 
sicrystal and another one with a lower mixing ratio. In 
the thermodynamic limit, we can neglect energy contri- 
butions from the interfacial boundaries between the two 
phases and optimize their respective lattice constants (or 
their 'volume fractions') to minimize the total energy per 
particle Ei of the two-phase state. 

Some plausible alternative structures are depicted in 
fig- 12 (''') two hexagonal lattices (phase separation of A 
and B dipoles), (ii) a. hexagonal A lattice with B dipoles 
in the triangle centers and a hexagonal A lattice, (in) a 
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FIG. 2: Alternative ordered structures and their energies per 
particle {Ei) compared to that of the quasicrystalline struc- 
ture Plotted is the fraction Ei/E^ as a function of the 
ratio of dipole moments D. 



hexagonal A lattice with B dipoles in every second tri- 
angle center and a hexagonal B lattice, (iv) a centered 
square lattice of A and B dipoles and a hexagonal lattice 
of B dipoles, and (v) a. centered square lattice of A and B 
dipoles together with the centered A-B hexagonal lattice 
of (ii). The corresponding values Ei for dipole strength 
ratios 1 < D < 10 are shown in the graph in fig. Irrelative 
to the value E^, of the quasicrystalline structure. 

The figure shows the preference of hexagonal order 
and the degeneracy of structures (i) and (ii) ioi D — 1, 
i.e., identical A and B dipoles. Furthermore, hexagonal 
ordering is preferred for _D < 3, whereas for larger D, 
partial tetragonal ordering seems energetically favorable. 
The quasicrystalline structure is closest to the optimum 
structure in the range 4 < Z? < 6, with an optimum 
around D ~ 4.5. However, the phase-separated struc- 
ture (v) clearly has lower energy. Accordingly, the qua- 
sicrystalline structure cannot be the ground state of the 
dipolar mixture. 

Still, we find that from a purely energetic point of 
view the quasicrystalline structure is a surprisingly com- 
petitive type of ordering for a binary mixture of dipoles 
within the appropriate mixing ratio and dipole strength 
ratio. In view of the fact that the alternative ordered 
structures (i)-(v) require a phase boundary, which pos- 
sesses a surface energy, the quasicrystal may become 
the preferred structure in finite systems. Interestingly, 
the quasicrystal becomes more favorable when the inter- 
action potential -Eint ~ 'r~^ is modified to a fictitious 
Eint ^ J' " with a close to 2 (a > 2 for reasons of con- 
vergence). Albeit we are not aware of any realization 
of such a potential in nature, this indicates that other 
types of scale- free interactions more strongly support the 
quasicrystalline ordering. 



IV. STABILITY ANALYSIS OF REFERENCE 
STRUCTURES 

A. Steepest descent calculations 

In this section, we test the dynamical stability of the 
quasicrystalline structure using the method of steepest 
descent. We displace the dipoles in small steps along 
their potential gradient with the step length being pro- 
portional to the modulus of the gradient vector. In nu- 
merics, this algorithm is used to find minima of multi- 
parametric functions |19| |. In physical terms, it describes 
the overdamped motion of the particles at zero temper- 
ature, which may be considered the simplest type of dy- 
namics to be implemented in a system. The method has 
been applied before to Lennard-Jones quasicrystals [2?il |. 

In our simulations, the maximum step length (for the 
dipole experiencing the largest force) is limited to a fixed 
value of order 0.01-0.001. Typically, after lO^-lO-* steps, 
the relaxation is finished when the dipoles start to oscil- 
late about fixed positions, where the amplitudes are of 
the order of the maximum step length. We perform two 
types of steepest descent calculations for varying dipole 
strength ratio D. In the first one, we start from the 
ideal quasicrystalline positions of the dipoles, and in the 
second one, we perturb the ideal positions by Gaussian 
random noise of different strength. 

The results of these simulations are somewhat hard to 
quantify since their concrete numerical outcome depends 
on the details of the algorithmic implementation, e.g., 
the value of the step length. The long-range dipolar 
potential and the lack of simple translational symmetry 
(which would let forces balance out trivially) give rise to 
a large number of dipole configurations being - mostly 
shallow - local energetic minima. Which of the different 
minima will be reached in a steepest descent calculation 
depends on computational details and initial conditions. 
Nonetheless, robust trends in the relaxation behavior are 
revealed and we obtain a reliable picture of the systems' 
dynamical stability. 

B. Stability as a function of dipole strength ratio D 

Figure O shows the order parameter as a function of 
the dipole strength ratio D as obtained after steepest de- 
scent calculations starting from the ideal quasicrystalline 
structure. The data indicate that the 'relaxed' config- 
urations stay closest to the ideal structure in a range 
4.5 < D < 5.5 with an optimum around D ~ 4.8. The 
optimum range of D essentially coincides with the range 
deduced from the comparison with alternative ordered 
structures in sec. IIIII 

This result is reproducible for different values of the 
step length and for different approximants (see fig. 
though the values of the relaxed order parameters </> may 
differ. This difference is not particularly significant, since 
the magnitude of (j) is very sensitive even to small dis- 
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FIG. 3: Order parameter <j> after steepest descent calcula- 
tions starting from the ideal quasicrystalline structure as a 
function of the dipole strength ratio D for three different ap- 
proximants. 

placements of the dipoles. For example, we find 4> — 0.8 if 
the dipole coordinates are randomly altered with respect 
to their values in the ideal quasicrystalline structure by 
a Gaussian noise with standard deviation a — 0.02 (cf. 
fig. 0. Therefore, it is reasonable to consider the 're- 
laxed' structures to represent the ideal quasicrystalline 
configuration. 

C. Recovery from perturbations 

Next, we test for metastability of the structure and ask 
if the quasicrystalline order recovers from small pertur- 
bations. Therefore, we apply Gaussian noise with zero 
mean and standard deviation a to each particle coordi- 
nate before letting them relax via steepest descent. The 
average order parameter (p of the final configurations is 
plotted in fig. 0] as a function of a for a few values of D 
from the optimum range. 

Up to an average initial displacement a ^ 0.15, the 
system relaxes back to the quasicrystalline structure, 
whereas for larger a, numerous defects remain as reflected 
by the decreasing average order parameter (f>. This be- 
havior is reminiscent of the empirical Lindemann crite- 
rion, which predicts the melting of a crystal for an aver- 
age displacement of 0.14 in units of the lattice constant. 
We note also that due to the initial Gaussian noise the 
order parameter decreases below 0.1 at a = 0.15 from 
which it recovers its optimum value around 0.9. This 
means that the quasicrystalline state is relatively robust 
against random perturbations. 

Figures 131 and 0] also include results for larger approx- 
imants. Clearly, the increasing system size does not lead 
to a significantly enhanced or reduced stability of the 
quasicrystalline order. For example, in Fig. |31 the order 
parameter at D = 4.8 is reduced to 0.85 in the medium 
approximant but recovers a value above 0.9 in the large 
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FIG. 4: Average order parameter (j) after steepest descent 
calculations starting from a perturbed structure for different 
values of D. (p is plotted as a function of the standard devi- 
ation a of the initial Gaussian displacements. As indicated, 
one curve corresponds to a larger approximant (1700 instead 
of 890 dipoles). For comparison, the order parameter (p due to 
pure Gaussian noise is shown, i.e., before the steepest descent 
relaxation starts. 

approximant. Hence we conclude that the stability of 
the quasicrystalline binary tiling has to be an intrinsic 
property of the ordered structure rather than a finite size 
effect. 



V. BEHAVIOR AT FINITE TEMPERATURES 

A. Monte Carlo simulations 

To assess the behavior of the binary system at finite 
temperatures we choose Monte Carlo (MC) simulations. 
While computationally less costly than e.g. Langevin dy- 
namics or even more detailed schemes, the dynamical MC 
simulations allow us to explore thermodynamical equilib- 
rium states and, with a reasonable choice of jump trials, 
should also yield a realistic scenario of the systems' evo- 
lution. 

As a jump trial in our simulations, a particle is chosen 
at random and imposed a Gaussian distributed displace- 
ment. The trial is accepted according to the Metropolis 
rule. The standard deviation of the Gaussian in the range 
0.001-0.1 is adjusted dynamically to ensure an efficient 
acceptance rate of trials in the range 10-60%. 

For sufficiently low temperatures, the outcome of these 
standard MC simulations goes well together with the re- 
sults from the steepest descent simulations. For example, 
starting from the ideal or slightly perturbed quasicrys- 
talline structure at T = 0.005, we find the order parame- 
ter first to decay and then to fluctuate around ~ 0.75, 
which confirms that (nearly) quasicrystalline order rep- 
resents a local energetic minimum in phase space. 

However, especially at slightly higher T, the degree 
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FIG. 5: Example trajectories of the order parameter (p as 
a function time (MC steps) from standard MC simulations 
at two different (fixed) temperatures T. The simulations are 
for D = 5 in an approximant of 760 dipoles. For each T, 
'run 1' started from the ideal quasicrystalline structure, while 
the initial configuration of 'run 2' was perturbed by Gaussian 
noise with t7 = 0.1 (cf. sec. Eg. 



and speed of the order parameter relaxation vary strongly 
from run to run in the MC simulations. In fig.|Sl example 
trajectories of are shown for T = 0.01 and 0.04. The 
quasicrystalline structure is relatively stable at T = 0.01, 
but the attained value of (f> depends erratically on the ini- 
tial conditions of the simulation. For example, in fig. [S] 
the order parameter emerging from an initially perturbed 
structure (run 2) unexpectedly exceeds the order param- 
eter emerging from the ideal quasicrystalline structure 
(run 1). Moreover, at T = 0.04 a slow and unpredictable 
decay of (p becomes observable on the time scale reached 
in the simulation. Even after excessively long runs of 
more than lO'^ MC steps, it remains unclear whether the 
system has reached an equilibrated state. 



method is to circumvent trapping of a systems' dynamics 
in local energetic minima at low temperatures by occa- 
sionally interchanging the configuration with the one of 
the same system simulated in parallel at higher temper- 
atures. Here, we consider 26 copies of our system at 
temperatures 2.5 x 10^"* = Ti > T2 > . . . > Tag = 0.071. 
Every 2000 MC steps, the particle configurations at ad- 
jacent temperatures T^, Ti_|_i are interchanged with a 
Metropolis type rate, 



exp 



{Ei+i — Ei) if Ei+i > Ei 
1 else . 

. 

where E^ and i^i+i are the energies of the configurations. 
It can be shown that this scheme allows different config- 
urations to occur with their correct Boltzmann weight at 
any of the temperatures Ti. 

In our implementation of the algorithm, a control pro- 
gram on a single PC keeps track of the configurations 
simulated at the various temperatures T^. Once the as- 
signment of a certain configuration to a Ti has been made 
for the next 2000 MC steps, it is passed as independent 
computing job to our queuing system. So the simula- 
tions can be carried out on a variable number of available 
CPUs which may also differ in speed. One MC step for a 
single configuration takes about 0.5 s on a contemporary 
Intel Pentium IV 2.8 GHz CPU. 

We find the parallel tempering algorithm to be highly 
effective in thermalizing our ensemble of 26 systems. This 
relatively large number allows us to keep the spacing be- 
tween the Ti small and thus to change configurations fre- 
quently while still covering a large range of temperatures 
from possible ordering to apparently fluid-like behavior. 
After typically 300 rounds (i.e. 300 x 2000 MC steps), 
we see no qualitative differences any more between the 
extremes of an ensemble started from the ideal quasicrys- 
talline structure and one started from random initial po- 
sitions. 



B. Parallel tempering 

To accelerate thermalization in the simulations, there 
are different possibilities. One of them is to introduce 
'artificial' multi particle flips as additional MC moves, as 
was previously done for the Lennard- Jones system in . 
In an extra series of our simulations, we tested an elemen- 
tary version of such flips, whereby every 100 MC steps 
the positions of an A and a B dipole are interchanged and 
the new configuration is evolved for several MC steps be- 
fore the whole move is either accepted or rejected. In 
general, we find the relaxation of </> to be faster, but the 
overall behavior remains unaltered. 

For the main part of our simulations, we use the stan- 
dard local moves of single particles, which might be closer 
to the dynamics of the experimental systems, but we ap- 
ply a more sophisticated thermalization scheme known 
as parallel tempering [21I |2^ . Basically, the idea of the 



C. Results: local ordering 

Figure|n|shows the mean order parameters (j) and 0„ for 
several n as a function of temperature T. The values are 
obtained from particle configurations at the respective Ti, 
each taken at the end of the 2000 MC step cycles of the 
parallel tempering scheme. The simulations are carried 
out at a dipole strength ratio D = 4.75 and started from 
random initial positions. The behavior is very similar for 
D — A.5 and 5. Apart from larger initial fluctuations 
of the order parameters (j) and 0io, the same holds for 
simulations started from ordered initial positions. 

In flg. \7\ part of an example configuration from the 
lower temperature range in the parallel-tempering en- 
semble is shown. In this representative example, no long- 
range quasicrystalline ordering is discernable. This find- 
ing is further supported by fig. El The order parameter 
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FIG. 6; Order parameters <j) and 0„ for several n as a function 
of temperature T as obtained from the parallel tempering MC 
simulations. The dipole strength ratio is D = 4.75. Note the 
jump in the ordinate. 



prevent the dipolar mixture from reaching its possible 
ground state as deduced for infinite system size. 

Based on these results, we can exclude the sponta- 
neous occurrence of long-range quasicrystalline order in 
a thermodynamically stable phase. Moreover, irrespec- 
tive of the (nearly) quasicrystalline structure being stable 
against small mechanical perturbations, we think that 
it may not represent a thermodynamically metastable 
state. Such a state would become separated from other 
states by an infinite free energy barrier in the limit of 
infinite system size. 

In accordance with experimental observations |[9|, [lOj , 
we find, however, that, while the overall structure is 
amorphous, there occur small domains with particles ar- 
ranged as in the quasicrystal, see fig. [3 As can be seen 
from the value of cf)^ and 0io compared to 0„ with n = 4, 8 
and 6, 12 in fig. El there is a tendency towards a preferred 
local 5- or 10-fold symmetry at lower temperatures. 
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FIG. 7: Configuration from the parallel tempering MC sim- 
ulations for D = 4.75. The temperature is T = 0.00134, the 
order parameter <f) ~ 0.07 and energy per particle E ~ 8.087. 
The part of the system shown is at the same scale as in fig.Q 
The lines highlight local quasicrystalline ordering and are 
drawn whenever the bond angles between nearest neighbors 
of A (B) particles are even (odd) multiples of 27r/10 (within 
an error of 10%). 



stays very small throughout the whole temperature range 
and there is no signature of a phase transition. 

We find that the mean potential energy per particle 
in the final configurations approximately fulfills E{T) ~ 
8.083 + T at low T, where the simple dependence on 
T can be understood from a harmonic approximation 
around the configurations with lowest energy. Thus, at 
low temperatures, typical configurational energies are be- 
low both the energy of the ideal quasicrystalline structure 
{Ei, ~ 8.103 for D — 4.75) and the one reached in the 
steepest descent calculations, cf. sec. IIVI On the other 
hand, E is still above the energy of the optimum phase 
separated structure (i?(i,) — 8.078) found in the ground 
state energy calculations in sec. IIIII This is not sur- 
prising, since the necessity of a phase boundary might 



VI. CONCLUSIONS 

Our results from the parallel tempering method pro- 
vide ample evidence that the quasicrystalline binary 
tiling with two sorts of dipoles does not correspond to 
a thermodynamical equilibrium state. Nevertheless, we 
can identify a range for the dipole strength ratio where 
a long-range quasicrystalline structure corresponds to a 
local minimum in the potential energy landscape of the 
system. This local minimum has an attraction basin cov- 
ering Gaussian fluctuations of the particle positions up 
to 15% of the distance of two neighboring A- and B- 
dipoles in the ideal reference structure. Our simulations, 
however, do not indicate that the barriers separating the 
local minimum from other minima increase with the sys- 
tem size. Therefore we do not expect that such a frozen 
state with quasicrystalline order exists although we can- 
not strictly exclude this possibility [23j . 

In any case, at low temperatures the times for the sys- 
tem to escape the local minimum by surmounting free 
energy barriers can become rather long. This is clearly 
seen in the kinetics modeled by conventional Monte Carlo 
simulations. Hence structures with long-range quasicrys- 
talline order can be kinetically stable over sufficiently 
long time to make them an interesting subject for fur- 
ther study and eventual applications. What "sufficiently 
long" in this context means could be tested in experi- 
ments. Todays optical tweezer techniques allow colloidal 
particles to be placed at defined positions. Accordingly, 
one could prepare a quasicystalline pattern and monitor 
its stability. Moreover, special boundary conditions and 
external stimuli may support the formation of quasicrys- 
talline structures. We also found that a modified scale- 
free interaction potential (x r^*^^+'^^ leads, for e > be- 
coming small, to ground state energies of the quasicrystal 
that within numerical error bars cannot be distinguished 
from the most favorable phase-separated lattice struc- 
tures investigated in section ITTll 
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In agreement with experiments we find that even in the 
disordered ground state of the binary dipole system, lo- 
cal bond orientational order with quasicrystalline 5- and 
10-fold symmetry is preferred. We made quantitative 
predictions for the temperature behavior of several local 
bond-orientation order parameters, which can be tested 
in experiments by using, for example, two-dimensional 
bi nary mixtures of superparamagnetic colloidal particles 

iliiliil- 



These symmetries are used to reduce the storage needs 
for the matrices of energy values and force components 
on the grid. 

For short notation, we use in this appendix Lj, = 1 as 
our length unit and ^QmiUij /AnLy as our energy unit, 
where rrii and nij are the magnetic moments of the two 
dipoles with pair vector Vij (transformation to the units 
used in the main text follows after elementary rescaling). 
Then we have 
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APPENDIX: ENERGY AND FORCE 
CALCULATIONS IN THE SIMULATIONS 

1. General procedure 

For an efficient calculation of the energy E{x,y) of (or 
force F(x, y) ~ — WE{x, y) acting on) a pair of dipoles i, j 
with distance vector r^j = {x,y), including their images 
in the periodically continued systems of the simulation 
box, we store the corresponding values on a fine grid 
of pair vectors and use a linear interpolation to obtain 
the values for r^j- in the continuum. According to the 
"minimum image convention" , possible distances in x- 
and y~ direction fall in the range —Lx/2 < x < Lx/2 
and —Ly/2 < y < Ly/2, where Lx and Ly arc the lengths 
of the system in the x~ and y-direction. We define by 
7 = Lx/ Ly the aspect ratio. Due to symmetry, E{x, y) is 
an even function of x and y, while the force components 
Fx{x,y) = -dxE{x,y), Fy{x,y) = -dyE{x,y) are odd 
functions of x, j/, and even functions of y, x, respectively. 



E{x,y)= 



1 



.=-oo [(x + 7Ai)2 + (y + !.)2 



3/2 



(A.l) 



The numerical calculation of this absolutely convergent 
series can be done by different means, for example by 
employing a two-dimensional variant of the Ewald sum- 
mation 24] or a simple extrapolation scheme, cf. We 
applied a method developed previously in our group , 
where the series in (|A.1|I is decomposed into an inner part 
for distances y^{x + 7/1)^ + (j/ + i^P smaller than a cutoff 
radius r,,!, and a remaining outer part, E = Ei^ + i?out- 
The inner part Ein is calculated by explicitly perform- 
ing the summation, while the outer part i?out is approx- 
imated by an integral. An analogous decomposition is 
done for the force components. In the following we dis- 
cuss the integral approximation for the outer parts and 
their numerical evaluation. 



2. Integral approximation for the energy 

Defining for a = 0, 1 



a else , 



where ceil(a:) is the lowest integer number larger than or 
equal to x, we obtain 



E„ 



E 

^——00 



(A.3) 



+\ ( [{^ + 7m)' + (2/ - + [{^ + li^f + {y + Mf^))'] 

dfJ.f{fi;x,y) . 



(A.4) 
(A.5) 
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Here [. . . ] stands for [(2; + + {y + J^)^] ■ The two integrals over v together yield 



/ 



dv 



M^^) - y 



l'l{^J') + y 



(A.6) 



Due to the piecewise definition of Va (n) , it is useful to split the remaining integral over /i in (|A.5p according to its 
boundaries, with 



dlJ^f{lJ-]x,y)= j dfif{fj,;x,y)+ f{^i;x,y) = Eout,i + Eout,2 

-'|/i|<rm/7 ■^lul>r,n/'y 



A'l>'"m/7 



The first part i?out,i is calculated analytically. With the abbreviations 



b±=\ 1 



y + 1 

r,„ ± X 



it reads 



En 



1 



2r„ 



2 — a_ — a I 2 — 6_ — & I 2 ^_ ^ 2 ^_ ^ 



y 



y + l 



2(y + i) 



If y = 0, the limit ?/ — > should be taken explicitly to avoid numerical instabilities, 



lim Sout.i = - 



1 



1 



1 / 1 



+ (2-6_ -6+) + 



1 



1 - 



1 1 



The second part -Eout,2 is approximated by a sum, 

Eont,2- X! ^{^^'^x,y) + e\[f{-^ll]x,y) + f{^ii]x,y)] , 



(A.7) 



(A., 



(A.9) 



(A.IO) 



(A.ll) 



where ^1 = ceil(r„i/7) and e = fix — ri-n/"f- For numerical stability, here the limit a; ^ of the addend with fi — 
should be considered explicitly. 



Hm /(O; x,y) = l- [{r„, - y) + (r,„ + y) ^ + (r„i - y) ^ + (r,„ + yY 



(A.12) 



In the numerics, the evaluation of (|A.lip can conve- 
niently be combined with the summation for the inner 
part Ein. For the energy matrix (and for the force calcu- 
lation discussed in the next section), we use = 10. If 
E is to be calculated repeatedly in a simulation, rm — 5 
might be a reasonable choice. In ample tests in the range 
0.8 < 7 < 1.2, we find the relative error of the integral 
approximation of the energy not to exceed 8 x 10"'' for 



Tm = 5 and 10 ^ for = 10. 

3. Integral approximation for the force 

The outer part of the x component force is 



y-Mf^) 



y + '^l{^^) 



2: + 7M \ [{x + -ft,Y + {y-Mt^Wf' [{x + jf^Y + {y + ^x{n)Yf\ 
I y-i^nin) y + vi{^.) \ 



{x + ^^if \[{x + ^^xY + {y-vo{^l)Yf'^ [{x + ^^iY + {y + u,{iYY]"^ ^ 

+ 7/i) / 1 ^ 1 ^ 

2 \[{x + ^iiY + {y-Mii)Yf^ [{x + ^^,f + {y + u^{p)Yf\ 

00 

dfJ.g{f^;x,y) . 



(A.13) 
(A.14) 
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Analogous to (|A.7p . the integral over /i is split into two parts, of which the first one yields 



out.l 



- - 

7 V 



1 



1 



- a- - 



y + l 



1 



1 



1 1 



- 6- - 



1 



a+ 



|?'m +x\ \b- 
1 /I 1 



1 1 ^ 



+2 

while the second part is approximated by a sum analogous to (|A.11|I . 

/ii-i 



^out,2 _ 9{^^■,x,y)+ (^-e\[g{-fii;x,y) + g{fii;x,y)] 



/i=-/ii+i 



(A.15) 



(A.16) 



Since the approximation turns out to be inaccurate in 
the neighborhood of x = and x = 0.57, it is advan- 
tageous to perform the force calculation in the ranges 
|x|/7 < 5 X 10-"* and 0.5 - |x|/7 < IQ-^ by explicit 
summation. 



The corresponding expressions for the force component 
can be obtained from the expressions for by 
interchanging x and y, replacing 7 by I/7, and rescaling 
by 1/74, i.e. F°-\x,y-^) - l/-f^ F°-\y/-f,x/r,l/l)- 
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